###reduced version
use.formula = as.formula(paste('ingroup.task ~ ',master.formula,sep = ''))
use.formula.out = as.formula(paste('outgroup.task ~ ',master.formula,sep = ''))
use.formula.reduced = as.formula(paste('ingroup.task ~ ',master.formula.reduced,sep = ''))

use.vars = c('ingroup.task','outgroup.task',master.use.vars)


#UO
use.dat.uo = dat.uo[,use.vars]
use.dat.uo = use.dat.uo[complete.cases(use.dat.uo),]

##N per city for weighted regression below
city.N = tapply(use.dat.uo$income,use.dat.uo$index.city,length)
city.N = as.data.frame(city.N)
city.N$index.city = row.names(city.N)
use.dat.uo = merge(use.dat.uo,
                   city.N)
use.dat.uo$index.city.numeric = as.numeric(use.dat.uo$index.city)


reg.uo.task.reduced = lm(use.formula.reduced, data = use.dat.uo)
clustered.se = cl(use.dat.uo, reg.uo.task.reduced, use.dat.uo$index.city)
reg.uo.task.reduced$se = clustered.se
reg.uo.task.reduced.weighted = lm(use.formula.reduced, data = use.dat.uo,
                                  weights = city.N)
brl.se <- brl(use.formula.reduced, dataframe=use.dat.uo, 
              clusterName="index.city.numeric", method="brl")[,2]
reg.uo.task.reduced.brl = reg.uo.task.reduced
reg.uo.task.reduced.brl$se = brl.se


reg.uo.task = lm(use.formula, data = use.dat.uo)
clustered.se = cl(use.dat.uo, reg.uo.task, use.dat.uo$index.city)
reg.uo.task$se = clustered.se
reg.uo.task.weighted = lm(use.formula, data = use.dat.uo,
                          weights = city.N)
brl.se <- brl(use.formula, dataframe=use.dat.uo, 
              clusterName="index.city.numeric", method="brl")[,2]
reg.uo.task.brl = reg.uo.task
reg.uo.task.brl$se = brl.se


##Non UO
use.dat.sec = dat.non.uo[,use.vars]
use.dat.sec = use.dat.sec[complete.cases(use.dat.sec),]

##N per city for weighted regression below
city.N = tapply(use.dat.sec$income,use.dat.sec$index.city,length)
city.N = as.data.frame(city.N)
city.N$index.city = row.names(city.N)
use.dat.sec = merge(use.dat.sec,
                    city.N)
use.dat.sec$index.city.numeric = as.numeric(use.dat.sec$index.city)

reg.sec.task.reduced = lm(use.formula.reduced, data = use.dat.sec)
clustered.se = cl(use.dat.sec, reg.sec.task.reduced, use.dat.sec$index.city)
reg.sec.task.reduced$se = clustered.se
reg.sec.task.reduced.weighted = lm(use.formula.reduced, data = use.dat.sec,
                                   weights = city.N)
brl.se <- brl(use.formula.reduced, dataframe=use.dat.sec, 
              clusterName="index.city.numeric", method="brl")[,2]
reg.sec.task.reduced.brl = reg.sec.task.reduced
reg.sec.task.reduced.brl$se = brl.se

reg.sec.task = lm(use.formula, data = use.dat.sec)
clustered.se = cl(use.dat.sec, reg.sec.task, use.dat.sec$index.city)
reg.sec.task$se = clustered.se
reg.sec.task.weighted = lm(use.formula, data = use.dat.sec,
                           weights = city.N)
brl.se <- brl(use.formula, dataframe=use.dat.sec, 
              clusterName="index.city.numeric", method="brl")[,2]
reg.sec.task.brl = reg.sec.task
reg.sec.task.brl$se = brl.se

###############
outtable = apsrtable(reg.uo.task.reduced,
                     reg.uo.task,
                     reg.sec.task.reduced,
                     reg.sec.task, 
                     Sweave = T,
                     coef.names = coef.names,
                     model.names = c('UO','UO','Non UO','Non UO'),
                     notes = ''
              )
writeLines(
  outtable, 'output/appendix/Table_A8.tex')


outtable = apsrtable(reg.uo.task.reduced.weighted,
                     reg.uo.task.weighted,
                     reg.sec.task.reduced.weighted,
                     reg.sec.task.weighted, 
                     Sweave = T,
                     coef.names = coef.names,
                     model.names = c('UO','UO','Non UO','Non UO'),
                     notes = ''
                )
writeLines(
  outtable, 'output/appendix/Table_A17.tex')

outtable = apsrtable(reg.uo.task.reduced.brl,
                     reg.uo.task.brl,
                     reg.sec.task.reduced.brl,
                     reg.sec.task.brl, 
                     Sweave = T,
                     coef.names = coef.names,
                     model.names = c('UO','UO','Non UO','Non UO'),
                     notes = ''
)
writeLines(
  outtable, 'output/appendix/Table_A20.tex')


###zelig for effects of segregation

##UO
reg.uo = zelig(use.formula, data = use.dat.uo,
		model = 'ls',
		cite = F)

reg.uo.out = zelig(use.formula.out, data = use.dat.uo,
               model = 'ls',
               cite = F)


#expected values
outputs.uo = matrix(nrow = length(dissim.scale.uo), ncol = 6)
for(i in dissim.scale.uo){
	row.i = which(dissim.scale.uo == i)
	outx = setx(reg.uo, outgroup.dissim.yeshiva = i,  outgroup.yeshiva = uo.homo.value,
		jerusalem = 0,
		demo1.sex = sex, 
		demo1.age = age, 
		demo1.ethnicity = ethnicity,
		demo2.left_right = ideology,
		income = income.value, 
		college = college.value, 
		nonjewish_pcnt = nonjewish_pcnt.value,
		immigrant = 0
	)
	sims = sim(reg.uo,outx, sims = n.sims)
	outs = sims$qi[[1]]
	outputs.uo[row.i,1] =  mean(outs)
	outputs.uo[row.i,2] = quantile(outs,low.conf)
	outputs.uo[row.i,3] = quantile(outs,high.conf)

	outx = setx(reg.uo.out, outgroup.dissim.yeshiva = i,  outgroup.yeshiva = uo.homo.value,
	            jerusalem = 0,
	            demo1.sex = sex, 
	            demo1.age = age, 
	            demo1.ethnicity = ethnicity,
	            demo2.left_right = ideology,
	            income = income.value, 
	            college = college.value, 
	            nonjewish_pcnt = nonjewish_pcnt.value,
	            immigrant = 0)
	sims = sim(reg.uo.out,outx, sims = n.sims)
	outs = sims$qi[[1]]
  outputs.uo[row.i,4] =  mean(outs)
	outputs.uo[row.i,5] = quantile(outs,low.conf)
	outputs.uo[row.i,6] = quantile(outs,high.conf)
	
  }



###Non UO
reg.sec = zelig(use.formula, data = use.dat.sec,
		model = 'ls',
		cite = F)
reg.sec.out = zelig(use.formula.out, data = use.dat.sec,
                model = 'ls',
                cite = F)

outputs.sec = matrix(nrow = length(dissim.scale.sec), ncol = 6)

for(i in dissim.scale.sec){
	row.i = which(dissim.scale.sec == i)
	outx = setx(reg.sec, outgroup.dissim.yeshiva = i, outgroup.yeshiva = sec.homo.value,
		jerusalem = 0,
		demo1.sex = sex, 
		demo1.age = age, 
		demo1.ethnicity = ethnicity,
		demo2.left_right = ideology,
		income = income.value, 
		college = college.value, 
		nonjewish_pcnt = nonjewish_pcnt.value,
		immigrant = 0
	)
	sims = sim(reg.sec,outx, sims = n.sims)
	out = summary(sims)$stats["Expected Values: E(Y|X)"]
	outs = sims$qi[[1]]
	outputs.sec[row.i,1] = mean(outs)
	outputs.sec[row.i,2] = quantile(outs,low.conf)
	outputs.sec[row.i,3] = quantile(outs,high.conf)

	outx = setx(reg.sec.out, outgroup.dissim.yeshiva = i, outgroup.yeshiva = sec.homo.value,
	            jerusalem = 0,
	            demo1.sex = sex, 
	            demo1.age = age, 
	            demo1.ethnicity = ethnicity,
	            demo2.left_right = ideology,
	            income = income.value, 
	            college = college.value, 
	            nonjewish_pcnt = nonjewish_pcnt.value,
	            immigrant = 0)
	sims = sim(reg.sec.out,outx, sims = n.sims)
	outs = sims$qi[[1]]
	  out = summary(sims)$stats["Expected Values: E(Y|X)"]
	outputs.sec[row.i,4] = mean(outs)
	outputs.sec[row.i,5] = quantile(outs,low.conf)
	outputs.sec[row.i,6] = quantile(outs,high.conf)
	}


	
###UO	
pdf('output/Figure_2_efficacy_seg_uo.pdf',
    width = pdf.width,
    height = pdf.height)    	
par(mar = c(4.5,3,.25,.25))
par(las = 1)

plot(0,0,
	ylim = ylims,
	xlim = xlims.dissim.uo,
	type = 'n',
	xlab = 'Segregation',
  ylab = '',
  cex.lab = cex.lab,
  cex.axis = cex.axis)

abline(h = 0,
		col = 'gray',
		lty = 2)
lines(dissim.scale.uo,
	outputs.uo[,1]- outputs.uo[,4],
	lwd = 4,
	col = 'black'
	)
lines(dissim.scale.uo,
	outputs.uo[,2] - outputs.uo[,4],
	lty = 3,
	lwd = 3,
	col = 'black'
	)
lines(dissim.scale.uo,
	outputs.uo[,3] - outputs.uo[,4],
	lty = 3,
	lwd = 3,
	col = 'black'
	)
	
dev.off()
	
##Non UO
pdf('output/Figure_2_efficacy_seg_non_uo.pdf',
    width = pdf.width,
    height = pdf.height)  	
par(mar = c(3,3,.25,.25))

par(las = 1)

plot(0,0,
	ylim = ylims,
	xlim = xlims.dissim.sec,
	type = 'n',
	cex.axis = cex.axis,
	cex.lab = cex.lab,
  xlab = '',
  ylab = '')

abline(h = 0,
		col = 'gray',
		lty = 2)
lines(dissim.scale.sec,
	outputs.sec[,1] - outputs.sec[,4],
	lwd = 4,
	col = 'black'
	)
lines(dissim.scale.sec,
	outputs.sec[,2] - outputs.sec[,4],
	lty = 3,
	lwd = 3,
	col = 'black'
	)
lines(dissim.scale.sec,
	outputs.sec[,3] - outputs.sec[,4],
	lty = 3,
	lwd = 3,
	col = 'black'
	)
	
dev.off()


#####################
###zelig for effects of homogeneity

##UO
outputs.uo = matrix(nrow = length(homo.scale.uo), ncol = 6)
for(i in homo.scale.uo){
	row.i = which(homo.scale.uo == i)
	outx = setx(reg.uo,  outgroup.yeshiva = i, outgroup.dissim.yeshiva = mean.seg.value,
		jerusalem = 0,
		demo1.sex = sex, 
		demo1.age = age, 
		demo1.ethnicity = ethnicity,
		demo2.left_right = ideology,
		income = income.value, 
		college = college.value, 
		nonjewish_pcnt = nonjewish_pcnt.value,
		immigrant = 0
	)
	sims = sim(reg.uo,outx, sims = n.sims)
	out = summary(sims)$stats["Expected Values: E(Y|X)"]
	outs = sims$qi[[1]]
	outputs.uo[row.i,1] =  mean(outs)
	outputs.uo[row.i,2] = quantile(outs,low.conf)
	outputs.uo[row.i,3] = quantile(outs,high.conf)

	outx = setx(reg.uo.out,  outgroup.yeshiva = i, outgroup.dissim.yeshiva = mean.seg.value,
	            jerusalem = 0,
	            demo1.sex = sex, 
	            demo1.age = age, 
	            demo1.ethnicity = ethnicity,
	            demo2.left_right = ideology,
	            income = income.value, 
	            college = college.value, 
	            nonjewish_pcnt = nonjewish_pcnt.value,
	            immigrant = 0)
	sims = sim(reg.uo.out,outx, sims = n.sims)
	out = summary(sims)$stats["Expected Values: E(Y|X)"]
	outs = sims$qi[[1]]
  outputs.uo[row.i,4] =  mean(outs)
	outputs.uo[row.i,5] = quantile(outs,low.conf)
	outputs.uo[row.i,6] = quantile(outs,high.conf)
	
  }


	
outputs.sec = matrix(nrow = length(homo.scale.sec), ncol = 6)
for(i in homo.scale.sec){
	row.i = which(homo.scale.sec == i)
	outx = setx(reg.sec, outgroup.yeshiva = i, outgroup.dissim.yeshiva = mean.seg.value,
	            jerusalem = 0,
	            demo1.sex = sex, 
	            demo1.age = age, 
	            demo1.ethnicity = ethnicity,
	            demo2.left_right = ideology,
	            income = income.value, 
	            college = college.value, 
	            nonjewish_pcnt = nonjewish_pcnt.value,
	            immigrant = 0
	)
	sims = sim(reg.sec,outx, sims = n.sims)
	out = summary(sims)$stats["Expected Values: E(Y|X)"]
	outs = sims$qi[[1]]
	outputs.sec[row.i,1] = mean(outs)
	outputs.sec[row.i,2] = quantile(outs,low.conf)
	outputs.sec[row.i,3] = quantile(outs,high.conf)

	outx = setx(reg.sec.out, outgroup.yeshiva = i, outgroup.dissim.yeshiva = mean.seg.value,
	            jerusalem = 0,
	            demo1.sex = sex, 
	            demo1.age = age, 
	            demo1.ethnicity = ethnicity,
	            demo2.left_right = ideology,
	            income = income.value, 
	            college = college.value, 
	            nonjewish_pcnt = nonjewish_pcnt.value,
	            immigrant = 0)
	sims = sim(reg.sec.out,outx, sims = n.sims)
	out = summary(sims)$stats["Expected Values: E(Y|X)"]
	outs = sims$qi[[1]]
	outputs.sec[row.i,4] = mean(outs)
	outputs.sec[row.i,5] = quantile(outs,low.conf)
	outputs.sec[row.i,6] = quantile(outs,high.conf)
	
  }

	
#UO	
pdf('output/Figure_2_efficacy_proportion_uo.pdf',
    width = pdf.width,
    height = pdf.height)	
par(mar = c(4.5,3,.25,.25))

par(las = 1)

plot(0,0,
	ylim = ylims,
	xlim = xlims.homo.uo,
	type = 'n',
	xlab = 'Outgroup Proportion',
  ylab = '',
	cex.axis = cex.axis,
  cex.lab = cex.lab
  )

abline(h = 0,
		col = 'gray',
		lty = 2)
lines(homo.scale.uo,
	outputs.uo[,1] - outputs.uo[,4],
	lwd = 4,
	col = 'black'
	)
lines(homo.scale.uo,
	outputs.uo[,2] -  outputs.uo[,4],
	lty = 3,
	lwd = 3,
	col = 'black'
	)
lines(homo.scale.uo,
	outputs.uo[,3] - outputs.uo[,4],
	lty = 3,
	lwd = 3,
	col = 'black'
	)

dev.off()
	
##sec
pdf('output/Figure_2_efficacy_proportion_non_uo.pdf',
    width = pdf.width,
    height = pdf.height)  	
par(mar = c(3,3,.25,.25))
par(las = 1)

plot(0,0,
	ylim = ylims,
	xlim = xlims.homo.sec,
	type = 'n',
  cex.axis = cex.axis,
	cex.lab = cex.lab,
  xlab = '',
  ylab = ''
  )

abline(h = 0,
		col = 'gray',
		lty = 2)


lines(homo.scale.sec,
	outputs.sec[,1] - outputs.sec[,4],
	lwd = 4,
	col = 'black'
	)
lines(homo.scale.sec,
	outputs.sec[,2] - outputs.sec[,4],
	lty = 3,
	lwd = 3,
	col = 'black'
	)
lines(homo.scale.sec,
	outputs.sec[,3] - outputs.sec[,4],
	lty = 3,
	lwd = 3,
	col = 'black'
	)
	
dev.off()

